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ABSTRACT 


An  analysis  and  computer  code  are  presented  to  calculate  the  effective 
wake  fraction  for  a  propeller  in  a  nonuniform  free  wake.  The  existence  of  an 
effective  wake  has  been  known  for  a  long  time,  but  the  rational  theory  pre¬ 
sented  here  seems  to  have  eluded  all  previous  investigators.  The  simple 
formula  developed  herein  to  calculate  the  effective  wake  requires  knowledge 
only  of  the  nominal  wake  and  the  load  distribution  on  the  propeller.  Since 
both  the  nominal  wake  and  the  load  distribution  are  ordinarily  given  in  the 
design  problem,  this  information  suffices  to  calculate  the  effective  wake 
for  a  given  design  once  and  for  all. 
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NOMENCLATURE 

A  constant  in  nominal  wake  distribution;  see  Eq.(22) 

B  constant  in  nominal  wake  distribution;  see  Eq • ( 22) 

CT  thrust  coefficient  based  on  ship  speed 
J  advance  ratio 

p  pressure 

AP  pressure  jump  on  actuator  disc 

r  radial  coordinate 

R  propeller  radius 

U(r)  nominal  wake  velocity 

U  effective  wake  velocity 

e 

u  axial  velocity 

u  induced  axial  velocity  at  actuator  disc 

o 

UQ  ship  speed 

V  reciprocal  of  U 

v  radial  velocity 

x  axial  coordinate 

F  pseudo-stream  function;  see  Eqs. (11) , (12) , (13) 

p  fluid  density 

cp  p/U 

stream  function  at  actuator  disc;  see  Eq,(2) 
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INTRODUCTION 

A  propeller  operating  in  the  wake  of  a  marine  vehicle  interacts 
with  the  wake  thereby  modifying  the  flow  into  the  propeller  from  that  of 
the  nominal  wake  measured  in  the  absence  of  the  propeller.  The  resulting 
change  in  the  velocity  distribution  at  the  propeller  plane  is  due  to  the 
change  in  pressure  induced  by  the  propeller  in  the  wake  flow  of  the 
vehicle.  The  radius  of  many  ship  propellers  is  about  the  same  order  of 
magnitude  as  the  thickness  of  the  wake  at  the  propeller  location.  The 
propeller  suction  therefore  increases  the  velocity  in  the  wake  in  the 
vicinity  of  the  propeller.  This  effective  inflow  velocity  distribution 
("effective  wake")  which  is  experienced  by  the  propeller  is  a  result  of 
the  interaction  of  the  propeller  and  the  wake  and  it  is  different  from  the 
nominal  velocity  distribution  (nominal  wake)  determined  in  the  absence  of 
the  propeller.  In  order  to  design  a  wake-adapted  propeller  for  specified 
propulsion  requirements  and  a  desired  radial  distribution  of  loading,  it 
is  important  to  have  an  estimate  of  the  effective  wake. 

1 

Huang,  et.al.,  have  shown  experimentally  that  the  conventional  tech¬ 
nique  of  adjusting  the  nominal  wake  by  a  constant  factor,  which  is  determined 
from  resistance,  self-propulsion,  and  open-water  experiments  such  that  the 
volumetric  average  effective  wake  is  equal  to  the  thrust  identity  wake,  is 
not  adequate  for  propellers  operating  near  the  stern  of  axisymmetric  bodies. 
Albeit  the  conventional  procedure  provides  a  reasonable  estimate  of  the 
average  inflow  to  the  propeller,  it  is  unable  to  determine  the  actual  radial 
distribution  of  this  inflow.  The  differences  between  the  effective  wake  and 
the  nominal  wake  were  found  (exper imentally)  to  be  greatest  near  the  propel¬ 
ler  hub  and  smallest  in  the  vicinity  of  the  propeller  tip.  This  result  is  in 
contrast  to  the  result  obtained  by  the  conventional  procedure,  which  usually 
yields  the  largest  differences  between  the  two  wakes  in  the  vicinity  of 
the  propeller  tip.  Therefore,  analytical  procedures  are  needed  to  provide 
a  more  accurate  description  of  the  interaction  between  a  propeller  and  its 

'Superior  numbers  in  text  matter  refer  to  similarly  numbered  references 
listed  at  the  end  of  this  report. 
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inflow.  This  is  especially  true  because  the  nature  of  the  effective 
inflow  into  the  propeller  disc  is  known  to  play  a  critical  role  in 
predicting  powering  and  cavitation  performance. 

The  invest igation  reported  by  Huang,  et.al.,  represents  the  most 
comprehensive  study  in  recent  years  on  the  steady  propeller-hull  inter¬ 
action  problem.  They  presented  a  prediction  method  which  obviates  the 
conventional  technique  of  determining  the  effective  wake  distribution 
into  the  propeller  operating  near  the  stern  of  an  axi symmetric  body. 

They  utilized  a  moderately-loaded  lifting-line  model  of  the  propeller 
to  calculate  the  propeller  induced  velocities  (the  model  was  originally 

developed  for  an  inflow  containing  no  vorticity;  see,  for  example, 

2 

Lerbs  ).  The  circumferential  average  velocity  distribution  is  taken  to 
be  given  two-diameters  upstream,  where  the  influence  of  the  propeller 
is  assumed  negligible.  The  flow  field  from  this  location  up  to  the 
propeller  disc  is  calculated  numerically  by  solving  a  simplified  form  of 
the  Euler  equation  of  motion  for  an  inviscid  fluid  containing  vorticity. 

The  propeller  induced  flow  field  displaces  the  streamlines  from  the 
original  position  when  the  propeller  was  not  present.  [t  is  the  re¬ 
distribution  of  the  upstream  vorticity  by  the  contraction  of  the  stream¬ 
lines  which  causes  an  additional  induced  velocity  field  and  it  is  this 
induced  velocity  which  is  the  effective  wake  fraction.  In  the  Huang, 
et.al.  procedure,  the  effective  wake  is  taken  as  the  inflow  into  the  pro¬ 
peller,  Since  the  effective  wake  fraction  is  not  known  a  priori, an  iterative 
scheme  was  required  in  which  the  propeller  influence  on  its  inflow  was  first 
calculated  as  if  there  were  no  effective  wake  fraction  at  the  disc,  and 
continued  until  the  prediction  of  effective  wake  converged.  The  comparison 
with  experiments  showed  good  agreement  with  the  lightly-  to  moderately- 
loaded  propellers  considered. 

Recently,  Goodman^  developed  a  linearized  axial  momentum  theory  for 
a  1 ightly-loaded  actuator  disc  operating  in  an  axially  directed  shear  flow. 
He  examined  the  effect  of  the  vorticity  contained  in  the  shear  flow  experi¬ 
enced  by  the  propeller  on  the  optimum  propeller  characteristics.  He  found 
that  the  velocities  induced  by  the  propeller  at  its  location  were  sub¬ 
stantially  altered  due  to  the  shear.  The  results  he  presented  were  for  an 
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optimum  disc  loading  only  although  he  presented  formulae  for  arbitrary 
disc  loadings. 

The  investigation  described  herein  has  the  purpose  of  developing  a 
theoretical  analysis  to  compute  the  effective  wake  which  is  assumed  to  be  a 
consequence  of  the  intensification  of  the  inflow  vorticity  by  the  extension 
of  the  vortex  tubes  caused  by  the  propeller  induced  flow.  In  other  words, 
a  theory  for  predicting  the  propeller  induced  velocity  field  is  developed, 
which  takes  into  account  the  vorticity  in  the  onset  flow.  Goodman*s  theory 
was  extended  to  compute  the  propeller  induced  velocity  field  due  to  a  pro¬ 
peller  with  an  arbitrary  radial  distribution  of  loading  operating  in  a 
shear  flow;  therefore,  all  the  assumptions  inherent  in  his  axial -momentum, 
lightly-loaded,  actuator-di sc  theory  apply  here  as  well.  The  induced  flows 
due  to  the  propeller  interaction  with  its  inflow  can  be  separated  into  two 
components:  one  due  to  the  disc  itself,  and  the  second  due  to  the  rearrange¬ 
ment  of  the  onset  flow  vorticity.  The  latter  is  the  effective  wake  frac¬ 
tion.  A  formula  for  the  effective  wake  on  the  disc  is  derived  for  an 
arbitrary  distribution  of  inflow  vorticity  which  may  be  evaluated  by  desk 
calculation.  This  formula  should  prove  useful  in  the  preliminary  design 
of  propellers  behind  axisymmetric  bodies  where  a  first  approximation  to  the 
effective  wake  is  required  quickly  (any  errors  incurred  by  the  lightly- 
loaded  actuator  disc  theory  assumption  are  expected  to  be  less  than  the 
other  uncertainties  in  preliminary  design  analysis  which  are  taken  into 
account  by  applying  error  margins  to  the  estimated  resistance  of  the  body, 
etc.).  However,  if  the  design  being  considered  is  heavily  loaded,  then  the 
prediction  by  the  formula  presented  herein  should  be  viewed  with  caution. 

The  designer  may  want  to  adjust  the  effective  velocity  distribution  computed 
herein  as  if  it  were  for  nominal  wake  until  its  volumetric  mean  matches  the 
thrust  identity  wake.  Then,  at  least  partially,  the  effect  of  shear  may  be 
accounted  for  in  a  quick  manner  as  is  often  desired  in  preliminary  design 
analysis. 

The  results  presented  show  that  without  shear  the  formulation,  and 

computed  results  reduce  to  the  linearized  actuator  disc  theory  of  Hough  and 

4 

Ordway  as  expected.  An  attempt  was  made  to  compare  the  results  of  the 
present  theory  with  the  measurements  presented  by  Huang, et .al . ,  but  a  large 
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discrepancy  resulted*  There  are  two  possible  reasons  for  this  discrepancy: 
The  theory  is  applied  at  the  propeller  disc,  whereas  the  measurements  were 
taken  somewhat  ahead  of  the  disc.  Near  the  disc  the  axial  velocity  is  in¬ 
creasing  rapidly  and  so  at  a  small  distance  ahead  of  the  disc  the  axial 
velocities  may  be  considerably  smaller  than  they  are  on  the  disc  itself. 

The  second  reason  for  the  discrepancy  is  that  the  measurements  (and  the 
theory)  of  Huang,  et.al.,  were  made  near  a  propeller  mounted  on  a  body 
of  revolution  as  in  a  submarine  or  torpedo.  The  body  with  its  incumbent 
boundary  condition  thus  contributes  to  the  modification  of  the  flow  near 
the  propeller.  The  present  theory,  on  the  other  hand,  is  for  a  propeller 
in  a  free  wake.  To  be  sure,  the  wake  will  have  been  generated  by  a  body, 
but  that  body  is  several  propeller  diameters  ahead  of  the  propeller  and 
therefore  does  not  play  a  role  in  the  interaction  of  the  propeller  and  wake. 
This  situation  is  more  likely  to  be  the  case  for  a  ship  rather  than  a  sub¬ 
marine  or  torpedo.  Thus,  what  is  required  to  test  the  present  theory  is  a 
set  of  experiments  specifically  designed  for  the  purpose. 

This  research  was  sponsored  by  the  Naval  Sea  Systems  Command,  General 
Hydromechanics  Research  Program,  and  administered  by  the  David  W.  Taylor 
Naval  Ship  Research  and  Development  Center  under  Contract  N00014-79-C-02 39, 
Davidson  Laboratory  Project  4679/056. 
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ANALYSIS 


The  computation  of  the  ’'effective  wake11  from  a  hydrodynami c  viewpoint 
involves  determining  the  modi f i ca t ions  to  the  flow  around  the  stern  of  a 
ship  into  the  propeller  due  to  its  operation.  On  an  axisymmetric  body  where 
the  propeller  is  in  the  wake,  the  disturbance  of  the  propeller  is  felt  no 
more  than  two-propeller  diameters  upstream  of1 the  propeller  planed  The 
required  input  data  for  the  determ \ nat ion  of  the  "effective  wake11  are  the 
nominal  wake  velocities.  These  data  can  either  be  measured  experimentally 
or  computed  analytically  for  a  particular  hull.  A  calculation  procedure 
for  the  flow  near  the  tail  region  of  a  body  of  revolution  has  recently  been 
presented  by  Geller.^  Axisymmetric  bodies  will  be  investigated  in  order  to 
focus  on  the  physical  nature  of  the  complex  interation  between  a  propeller 
and  a  thick  wake.  Their  geometric  simplicity  offers  considerable  computa¬ 
tional  convenience  in  treating  the  fundamental  aspects  of  the  interaction 
physics. 

In  this  investigation  the  nominal  wake  is  assumed  to  be  given.  In 
addition,  the  nominal  wake  is  usually  supplied  in  the  design  problem  of 
the  propeller.  The  velocity  field  induced  by  the  propeller  is  to  be  calcu¬ 
lated  by  using  the  theoretical  procedure  described  next. 


Induced  Axial  Velocities  in  the  Plane  of  the  Propeller 

It  is  shown  by  Goodman  that  the  pressure  jump  across  the  propeller 
disc  in  a  nonuniform  stream  U(r)  is  given  by 


AP 

2p 


Uu  +  ~  >L 
o  r  dr  o 


where 


0) 

(2) 


defines  a  Stokes  stream  function,  and  is  the  perturbation  axial  velocity 
in  the  plane  of  the  propeller.  Equations  (1)  and  (2)  give  a  relationship 
between  the  pressure  jump  AP  and  the  induced  velocity  uq.  Our  object  will 
be  to  invert  the  relationship  so  that  the  induced  velocity  is  expressed 
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directly  in  terms  of  the  pressure  jump.  Solving  for  i|t  ,  we  obtain 


Ap 


t  5  -  [  ru  dr  =  ~ 

o  J  o  1 


-  u 


_  J_  du 

U  r  dr 


(3) 


Let  |  U  3r  "  ^(r)  9  TpIT  =  9  t*1en  upon  fferent i a t i ng  Eq.(3)  with 

respect  to  r  ,  there  is  obtained 


du 
_ _ c 

dr 


[f1  +  rf  ju  =  fg'  -  gf 1 


This  is  a  first  order  linear  differential  equation  for  u^,  and  its  solution 
can  be  shown  to  be 

AP  1  dU  pr  APrdr 


Uo  2pU  + 


au  p 
r  dr  J 


o  2P  IT 

Thus  the  total  velocity  in  the  plane  of  the  propeller  is 


(5) 


U  +  u  =  u  +  AL.  +  1  ^  f  ££#■«•, 
o  2pU  r  dr  J  2p  US 


(6) 


The  stream  velocity  U  is  frequently  called  the  nominal  wake,  while  the 
second  term  on  the  right-hand  side  is  the  induced  velocity  in  the  absence 
of  shear.  The  third  term  is  the  ihduced  velocity  due  to  the  presence  of 
shear,  and  it  is  customary  in  propeller  design  to  lump  the  first  and  third 
terms  together  and  call  them  an  effective  wake.  Thus, 

r  APrdr 


u  =u+t^.r 
e  r  dr  i  2p  u2 


(7) 


Note  that  in  a  uniform  flow  the  effective  and  nominal  wakes  become  identical. 
From  Eq.(7)  it  is  possible  to  calculate  the  effective  wake  directly  from  the 
load  distribution  on  the  propeller.  In  the  design  problem  the  load  distri¬ 
bution  is  ordinarily  given  so  that  the  effective  wake  can  be  calculated  for 
a  given  design  once  and  for  all.  This  may  be  contrasted  with  the  procedure 
of  Huang,  et.al.,^  in  which  the  effective  wake  is  calculated  from  flow 
variables  which  are  unknown  at  the  outset  so  that  it  becomes  necessary  to 
iterate  between  the  effective  wake  and  the  wake-adapted  propeller  to  which 
it  is  designed.  Equation  (7)  is,  in  fact,  so  simple  that  it  is  easy  to 
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calculate  the  effective  wake  by  hand  using  a  trapezoid  rule  to  carry  out 
the  integration.  The  simplicity  inherent  in  Eq.(7)  comes  about  because 
the  equations  of  motion  have  been  linearized  in  developing  Eq.(l),  and 
also  because  the  nominal  wake  has  been  taken  to  be  independent  of  the  axial 
coordinate.  These  assumptions  make  Eq.(l)  so  simple  that  it  can  be  in¬ 
verted*  By  contrast,  Huang,  et.al.,  have  not  linearized  and,  consequently, 
no  such  inversion  is  possible  in  their  case.  On  the  other  hand,  the  theory 
on  which  the  propeller  design  is  based  is  invariably  a  linear  one  so  that 
no  great  advantage  is  obtained  in  retaining  nonl i near i t ies  in  the  determina¬ 
tion  of  the  effective  wake.  With  regard  to  the  second  assumption,  viz., 
that  the  nominal  wake  is  independent  of  the  axial  coordinate,  it  is 
certainly  true  that  the  nominal  wake  varies  axially  according  to  boundary 
layer  theory  and,  if  the  Reynolds  number  is  sufficiently  large,  the  varia¬ 
tion  will  be  very  small  over  a  distance  of  one  or  two  propeller  diameters. 
Thus  it  is  sufficiently  accurate  to  take  the  nominal  wake  to  be  independent 
of  the  axial  coordinate  and  equal  to  its  value  in  the  plane  of  the  propeller 
in  the  absence  of  the  propeller.  With  these  considerations  it  can  be  seen 
that  the  approximations  inherent  in  the  development  of  Eq.(7)  are  no  more 
severe  than  the  approximations  at  the  core  of  the  method  of  Huang,  et.al., 
although  they  are  different. 

Of  course, the  present  theory  is  valid  for  a  propeller  in  a  free  wake 
(as  for  a  ship),  whereas  the  theory  of  Huang,  et.al.,  is  valid  for  a 
propeller  mounted  on  a  body  of  revolution  (as  for  a  submarine  or  torpedo). 
Thus  it  is  not  possible  to  evaluate  the  present  theory  using  the  data  of 
Huang,  et.al,  and  what  is  needed  is  an  experiment  specifically  designed 
to  test  the  present  theory. 


Induced  Axial  Velocities  Ahead  of  the  Propeller 

In  order  to  be  able  to  calculate  the  induced  velocities  ahead  of  the 
propeller  disc,  it  is  necessary  to  consider  the  complete  linearized  field 
equations.  The  linearized  Euler  equations  ahead  of  the  disc  in  the  case  of 
axial  symmetry  are 

Su  .  dU  _  IS 


..  Ou  ^  dU  _  1  dp 

Ox  v  dr  p  3x 


(8) 


37 
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7  3r 


(9) 
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The  equation  of  continuity  is 

(™)  +  (rv)  -  0  (10) 


In  order  to  develop  an  equation  for  the  axial  perturbation  velocity 
in  the  field,  we  define  a  pseudo-stream  function  T  such  that 


2  =  u 1  dr 

p  ~  r  dr 

do 

“--■?§?<ur> 

(12) 

v  ■  t  It  (ur) 

(13) 

With  this  definition,  Eqs.(8)  and  (10)  are  satisfied 
Eq. ( 12) ,  we  see  that 

identical ly. 

Expandi ng 

u  ar  r  du 
u  r  dr  "  7  dr 

=  _  1  £  _  £  dU 

U  p  r  dr 

0*0 

or,  from  Eq. (1 1) 

,  _  ,  ...  r  £  rdr 

_  1  £  _  l  dU  r  _ 

Up  rdr  j  yf 

(15) 

Equation  (15)  may  be  looked  upon  as  a  generalization  of  Eq.(5)  to  points 
ahead  of  the  propeller.  Indeed,  on  the  propeller  disc  p  =  -  j  AP  ,  and, 
upon  substituting  into  Eq .(15),  it  can  be  seen  that  this  equation  reduces 
identically  to  Eq.(5).  From  (14)  or  (15),  the  effective  wake  can  be  seen 
to  take  one  or  the  other  of  the  following  two  forms 

«  -u.i «/  K- 

e  r  dr  J  \j2 

o 

(16) 

-u-i&r 

r  dr 

(17) 

The  effective  wake  can  therefore  be  determined  ahead  of  the  propeller  once 
p/p  or  T  are  known  ahead  of  the  propeller.  One  possibility  is  to  use 

pressure  measurements  ahead  of  the  propeller  and  to  substitute  such 
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measurements  into  Eq.(l6).  We  will  be  concerned  here,  however,  with  an 
analytical  procedure,  and  for  this  purpose  will  require  the  solution  of 
the  field  equation  for  either  T  or  p/p.  First  consider  the  field  equation 
for  r . 

Upon  substituting  Eqs.(ll)  and  (13)  into  Eq.(9),  there  is  obtained 


In  order  to  obtain  the  field  equation  for  p  ,  we  differentiate 
Eq.(IO)  with  respect  to  x  and  substitute  Eqs.(8)  and  (9).  This  yields 


f_P/£  + 

dx2 


dU  dv  1  dU  dp/p 
dr  3x  ~  IT  dr  dr  ~  0 


(18) 


Then  using  Eq.(9)  again  to  eliminate  the  following  field  equation  for  p 
i s  found 


'fp  -  2 


uj 

u 


dP 

37 


r-  ¥•  =  o 


(19) 


This  is  basically  the  same  equation  for  the  pressure  found  by  von  Karman  and 
•  o 

Tsien.  It  can  be  manipulated  into  another  form  more  suitable  for  our  pur¬ 
poses.  Let 


(20) 


It  can  then  be  shown  that  Eq .(19)  reduces  identically  to  the  following 
equation  for  cp ; 

V  V^cp  -  epV5  V  =  0  (21) 

and  this  equation  for  cp  can  be  shown  to  be  valid  whether  or  not  U  and  cp 
also  depend  on  the  azimuthal  angle  as  well  as  the  radial  coordinate. 

We  will  choose  to  solve  Eq • (21)  rather  than  either  Eqs.(17)  or  (19) 
because  for  a  particular  family  of  wakes  it  admits  of  a  relatively  simple 
solution.  It  can  be  seen  from  Eq .(21)  that  whenever  the  wake  function  V 
obeys  the  Laplace  equation,  the  pressure  function  cpwill  also  obey  the 
Laplace  equation.  In  the  case  of  axial  symmetry  the  solution  of  the  Laplace 
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equation  that  V  must  satisfy  is  of  the  form 


1  =  V  =  A  4n  r  +  B 


(22) 


where  A  and  B  are  constants.  It  will  be  assumed  for  the  purposes  of  this 
development  that  the  nominal  wake  is,  indeed,  of  the  form  of  Eq.(22).  In 
this  case,  since  cp  obeys  the  Laplace  equation,  we  can  immediately  write 
down  the  solution  for  the  pressure  in  the  field  from  Lamb^  (Article  102, 

n) : 

P  =  -  I  I  J0(kr) J0(ks)e"klX*  Ay ['l}  ksdsdk  (23) 


where  AP  is  the  jump  in  pressure  across  the  disc  (=  AP(s)  for  s  <  b, 
*  0  for  s  >  b) . 

Thus  from  Eq.(l6),  the  effective  wake  is 


ue  =  u  +  77^7  I  dri  I  ds  I  dk 


ft.P.(s). 


Now  the  r^  integral,  viz., 


sr,kJo(krj)Jo(ks)e 

U(s)u(r) 


kx 


(24) 


r  Jo(kri)ridri 

“  O  «(r,) 

can  be  evaluated  explicitly  when  U  obeys  Eq.(22),  and  the  result  is  (see 
Appendix  A) : 


I 


r  Jj(kr)  A[1  -  JQ(kr) ] 


k  U(r)  k2 

Thus  the  effective  wake  becomes 


(25) 


Ue  =  U 


I  AU(r)  j? 

"  2  r  i  puTsl 


.  1  AV(r)  !?  ^(s)sl2(s,r)ds 

+  2  — r  - pu(i7 - 


(26) 
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where 


I,  ■  J  J0(ks)jj(kr)ekxdk 


*2  =  I 


"  Jo<kr)) 


JQ(ks)ekxdk 


and  use  has  been  made  of  the  relationship 

which  follows  from  Eq.(22). 

It  can  be  shown  by  using  identities  involving  the  Bessel  functions 
j  r 

that  *2  =  7  I  anc*  t^)en>  passing  to  the  limit  x  -*  0 ,  i.e.,  on  the  disc, 

°  o 

from  Gradshteyn  and  Ryztik  (#6*5123): 


•l"7  r>s 

=  0  r  <  s 


x  =  0 


and,  consequently, 


I  0 - j£n  —  r  >  s 

2  r  s  ^  a  9 


x  =  0 


1^  may  also  be  written 


,  1  (  1  1  ^ 

2  =  At  VUX7T  '  T3TsT^  r  >  s 


x  =  0 


r  <  s 


and  in  this  form  it  is  easy  to  show  that  Eq.(26)  reduces  identically  to 
Eq.(7)  on  the  disc  as  it  must. 
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COMPUTER  CODE 

A  FORTRAN  program  was  developed  to  compute  the  velocities  given  by 
Eqs.(7)  and  (26)  depending  upon  the  location  of  the  field  point  of  interest. 

If  the  field  point  is  on  the  disc,  i.e.,  at  x=0  and  any  r  ,  then  Eq.(7) 
is  used  and  U=U(r)  may  be  any  arbitrary  distribution  of  inflow  velocity 
(nominal  wake).  However,  if  the  field  point  is  off  the  disc,  i.e.,  x/0  and 
any  r,  then  the  inflow  velocity  distribution  must  be  approximated  by  Eq.(22). 
The  input  radial  distribution  of  the  thrust  must  gradually  drop  to  zero  at 
the  hub  and  tip  of  the  disc,  i.e.,  a  disc  with  constant  load  may  not  be  con¬ 
sidered.  This  is  not  a  severe  restriction  because  the  latter  is  unrealistic. 
This  restriction  is  implicit  in  the  assumptions  made  to  develop  Eq .  (26) ;  see 
Lamb.  ^ 

Appendix  B  presents  a  listing  of  the  code.  The  input/output  informa¬ 
tion  is  described  next. 

INPUT 

The  first  input  card  is  read  in  the  main  program  VFSHEAR;  see  page  B1« 
The  constant  to  be  read  is  an  integer  called  NCASES  with  format  15.  It  in¬ 
dicates  the  number  of  sets  of  input  cards  to  follow  this  card,  i.e.,  it  is 
the  number  of  cases.  One  case  constitutes  the  complete  set  of  input  cards 
read  in  the  subroutine  INPUT;  see  page  B3  in  Appendix  B. 

The  set  of  INPUT  cards  for  each  case  is  clearly  presented  in  the 
FORTRAN  statements  on  page  B3.  The  appropriate  formats  are  given  as  well. 
Therefore,  the  only  additional  l n formation  required  to  execute  the  program 
is  the  definition  of  each  piece  of  input  data.  The  FORTRAN  variables 
representing  the  input  data  are  defined  in  Table  1. 
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TABLE  1 


INPUT  TO  VFSHEAR 

DEFINITION  OF  FORTRAN  VARIABLES  IN  THE  ORDER  THEY  ARE  READ 


FORTRAN  VARIABLE 

TYPE 

DESCRIPTION 

1  DENT ( 1 ) 

Alphanumeric 

Any  phrase  used  to  define  the  case  to  be  calculated. 

RHOF 

Real 

Fluid  densi ty ,  p 

DS 

Real 

Annular  strip  width  (taken  to  be  0.05  if  the  hub 
radius  is  0.2R) . 

AWK 

Real 

The  wake  parameter  A  in  Equation  (22). 

CTHRUST 

Real 

The  thrust  coefficient,  C ^ 

UNORM 

Real 

Volumetric  average  of  the  input  velocities  (use 

1.0  if  field  points  are  off  the  disc) 

NSTRIP 

1 nteger 

Number  of  annular  strips  (16). 

s(D 

Real 

Radial  coordinate  of  the  center  of  each  strip. 

US(I) 

Real 

Axial  inflow  velocity  at  each  strip. 

DELPS(I) 

Real 

Pressure  jump  across  each  strip. 

NOFPS 

Integer 

Number  of  field  points  to  be  calculated. 

/  i 

N 

I nteger 

Number  of  Laguerre  coefficients  in  Guass-Laguerre 
quadruture  scheme  (50). 

XF(J),RF(J) 

II 

Rea  1 

Coordinates  of  the  field  points  of  interest. 

UR(  J) 

Real 

Axial  inflow  velocity  at  the  field  point. 

OUTPUT 


First,  the  number  of  cases  is  printed  followed  by  a  review  of  the  input 
data  for  the  case  under  cons i derat  ion.  UR  is  read  as  an  integer  if  the  point 
is  on  the  disc.  In  addition,  when  computing  points  on  the  disc,  the  field 
points  must  be  selected  such  that  they  are  at  the  center  of  the  strip.  The 
integer  value  of  UR  is  the  number  of  the  radial  strip  at  which  the  induced 
velocity  is  to  be  computed.  UR  ®  1.0  corresponds  to  the  strip  closest  to 
the  hub,  while  UR  =  16.0  corresponds  to  the  strip  nearest  the  tip.  If  the 
field  points  are  off  the  disc,  then  the  actual  velocity  is  input  for  UR. 

The  first  set  of  computed  output  is  called  ''velocity  at  given  field 
points"  and  is  calculated  using  Eq.(26).  If  AWK  =  0.0,  then  the  values  of  U 
do  not  account  for  shear  and  are  the  potential  part  of  the  propeller  induced 
flow. 

The  second  set  of  output  apply  only  when  the  field  points  are  taken  on 
the  disc.  The  formula  used  to  calculate  the  second  set  of  U  is  Eq.(7). 

If  the  program  is  executed  for  an  arbitrary  distribution  of  U(r)  for 
field  points  on  the  disc,  AWK  should  be  set  to  zero.  Then,  if  the  first 
set  of  U^  are  substracted  from  the  second  set  of  U's,  the  user  will  obtain 
the  linearized  actuator  disc  theory  prediction  of  the  effective  wake  frac¬ 
tion. 

For  field  points  off  the  disc,  the  effect  wake  fraction  is  determined 
by  running  two  cases.  In  the  first  case  AWK  is  set  to  its  appropriate 
value  for  the  input  US ( I ) .  In  the  second  case  everything  else  is  kept  the 
same  but  AWK  is  set  to  zero.  The  differences  in  the  finite  values  of  the 
output  U's  is  the  effective  wake  fraction. 
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RESULTS 


Since  the  present  theory  is  for  a  lightly-loaded  actuator  disc, 

for  U  =  constant,  we  should  recover  the  results  calculated  and  reported 

4 

by  Hough  and  Ordway.  The  results  calculated  with  the  present  code 
using  a  loading  distribution  close  to  Hough  and  Ordway*s  variable  distri¬ 
bution  were  compared  with  the  data  presented  by  Hough  and  Ordway  and  were 
found  to  be  virtually  identical.  Differences,  which  are  in  the  third 
significant  figure,  are  due  to  the  numerical  integration  scheme  and  also 
because  the  load  distribution  that  was  used  accounts  for  a  hub  while 
Hough  and  Ordway's  distribution  does  not. 

Next,  the  nominal  wake  presented  in  Figure  I6e  of  Huang,  et.al., 
was  used  together  with  the  appropriate  values  of  C^s(*a0.371)  and  J(=1.25) 
and  the  load  distribution  of  Hough  and  Ordway  (modified  to  account  for  a 
hub)  to  calculate  the  total  axial  velocity  using  Equation  (5)#  It  is  to 
be  noted  that  the  nominal  wake  was  measured  one-quarter  of  a  radius  ahead 
of  the  propeller  disc,  but  for  the  purposes  of  this  calculation  it  was 
assumed  that  the  nominal  wake  remained  the  same  at  the  disc  itself.  No 
attempt  was  made  to  use  the  off-disc  computer  program  for  this  comparison 
because  the  wake  measured  and  presented  by  Huang,  et.al.  does  not  fit 
Eq .  (22) .  Instead,  the  wake  appears  to  fit  the  "law  of  the  wall,"  viz., 
Usainr+b,  and  this  fact  was  used  to  extrapolate  the  nominal  wake  data 
presented  to  the  hub  of  the  propeller.  The  results  are  presented  in  Fig¬ 
ure  1.  A  set  of  experiments  specifically  designed  to  test  the  present 
theory  is  clearly  required.  Since  it  is  impossible  to  measure  the  axial 
velocity  directly  at  the  propeller  disc,  it  would  be  necessary  to  measure 
it  a  small  distance  ahead  and  an  equal  distance  behind  the  disc.  Then, 
since  the  axial  velocity  is  known  to  be  monotonic  and  continuous  through 
the  disc,  its  value  at  the  disc  can  be  inferred  by  averaging  the  fore  and 
aft  measurements.  Moreover,  in  order  to  simulate  the  wake  of  a  body,  a 
wake  screen  could  be  used  and  the  measurements  taken  in  a  water  tunnel 
using  a  Lazer-Doppler  Anemometer. 
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NOMINAL  WAKE  AND  TOTAL  INFLOW  VELOCITY  FROM  FIGURE  l6e 
OF  REFERENCE  C 1 3  AND  EQUATION  (6) 
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APPENDIX  A 

EVALUATION  OF  A  CERTAIN  INTEGRAL 

Consider  the  integral 

.  rr  Jo(krl>rldrl 

'4 

where 

^  =  A  in  Tj+B 
=  A  in  krj+B1 

where 

B'  =  B-A  in  k  . 

Let 

krl  =  Z1 
kr  -  z 

Then 

t  2 

1  ■  Ts  J  z^ZjCA  In  z,  +  B']Jo(z]) 
k  o 

but 

zlJo(2l>  •  • 

Upon  substituting  for  J^,  the  term  proportional  to  B  can  be  integrated 
directly,  while  the  term  proportional  to  A  can  be  integrated  by  parts. 
This  yields 

1  *  ~2  {azJ  j  (z)  -in  z-A  J  (Zj)dZj  +  b'zJj(z)} 

But  since  the  integral  can  be  evaluated.  Then,  upon  substituting 

for  B1  in  terms  of  B  and  simplifying,  the  result,  Equation  (25),  follows. 
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APPENDIX  B.  FORTRAN  LISTING  OF  PROGRAM  VFSHEAR 

PROGRAM  VFSHE  AR(  IN  PUT* OUTPUT,  TAP  b  5 -INP  UT*  T  A  P  E6-0UT  PUT*  T  OP  6  7* 
*  TAPE1,TAPE3) 

C 

C  PROGRAM  TO  CGMPUTE  VELOCITY  FIELD  UPSTREAM  OF  A  PROPELLER 

C  OPERATING  IN  A  SHEAR  FLOW  DEFINE  BY  I/O  «  AWK*LOGE(R)  ♦  BWK 

C 

COMMON/ INPUT /R HOF* S< 32 )* DS, XF < 10U ) , R F ( 100 ) , DE L P S ( 32 ) 
C0MM0N/NCMWAKE/US<32),NCFPS*UR(100),U<100), AWK 
C0MM0N/INTGRND/FI0,FU,FI2,XI0,XI1»  XI2 
DIMENSION  FI0(960)*F11(960)*F  12(960) 

DIMENSION  XI0(32)*XI1(32)*XI2(32) 

C0MM0N/LAG2RS/  Y(15)*N- 
COMMON/DUDR/  URX(IOC) 

COMMON/THRUST/  C T HROS T» UNOR M 
COMMON/STRIP/  NSTRIP 
COMM ON /LAG/ XL* AL* NN 
DIMENSION  Q(125), 01(125) 

DIMENSION  XL(50),AL(50)»B(50),C(50) 

C 

C  READ  INPUT  DATA 

C 

RE AD ( 5*  2G00 )  NCASES 
2000  FORMAT (IS ) 

WR I T E ( 6*  2002  )  NCASES 

2002  FORMAT  (  /  /  /*  NUMBER  OF  CASES  «*15/) 

DO  2001  INC I»1*NC ASES 
CALL  INPUT 

C  FOR  EACH  FIELD  POINT  COMPUTE  VELOCITY 

C** 

IF(USd).Ea.O.O)  8WK»US  (  2  ) 

I F  <  US ( 1 ) »NE • 0  •  0 )  GO  TO  95 
DO  96  IXN»1, NSTRIP 

96  US(IXN)«1.0/(AWK*AL0G(S(1XN)  )+BwK) 

95  CONTINUE 
C 

C  CALCULATION  OF  ZEROS  AND  COEFFICIENTS  IF  N . NE . 1 3 . OR . N. N E . 10 . 

C 

EPS  ■  0.0000000001 
EPS  •  EPS**2 

IF(N.EQ.15.0R.N.EQ.10 )  GO  TO  546 
NN  «  N 

DO  547  I  -  1 *  NN 
XI  -  FLOAT ( I ) 

B  (  I )  -  2.*XI-1. 

547  C  1 1  )  »  C  XI— 1 •  )**2 

CALL  LAGL'ER{NN,XL*AL*0.0*3*C*EPS*CSX*CSA*TSX*TSA) 

546  CONTINUE 

DO  1001  I FP* 1 »  NQF  PS 
IF ( XF  t IFPJ.EG.O.O)  GO  TO  1019 
C 

C  CALCULATE  INGREDIENTS  TO  VELOCITY  INTEGRATIONS 

C 

IFPX  «  IFP 

CALL  INTGRD(IFPX) 

C 

C  PEPFORM  GAUSS  I AN-LAGUERRE  QUADRATURE 

C 

CALL  GAUSSLG 
C 

C  AT  THIS  POINT  WE  HAVE  FOR  EACH  STRIP  10,  II, AND  12 

IF(XF(IFP).GT.0.2)  GO  TO  2090 
DO  10001  1X«1, NSTRIP  (Bl) 
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ARGQ*(XF(IFP)+*2+RF(IFP)**2+S(lX)**2)/(2.*S(IX)*RF(IFP)  ) 
FACQ»XF(1FP)/((RF(IFP)*S(IX))**1.5) 

CALL  LEG2HG(ARGQ»2.Q>Q1.2) 

XIO< IX )«Q1(1 )*FACQ/3. 1415927 
XIO(  IX)— XIO(  IX) 

10001  CONTINUE 
2090  CONTINUE 
C 

C  CALCULATE  VELOCITY  AT  PARTICULAR  FIELD  POINT 

C 

C  SPANWISE  INTEGRATION 


USUM1-0.0 

USUM2  •  0.0 
DO  1002  JS-l/NSTRIP 

FAC  •  DS*DELPS< JS > / C R HOF +US ( J S ) *2 . ) 

UINTG  »  SCJS)*FAC*XIO(JS) 

USUMO  ■  USUMO  +  UINTG 
UINTG1  -  S( JS)*FAC*XI1( JS) 

USUM1  •  USUM1  ♦  UINTG1 
UINTG2  «  S(JS)*FAC*XI2(JS) 

1002  USUM2  •  USUM2  ♦  UINTG2 

U( IFP)  ■  USUMO-AWK+UR  ( I F  P )  +  US  JM  /  P  F  ( I F  P  )  ♦  AWK*A  WK  *UR  (  IF  P  )  *  UR  (  I F  P  )  * 
1  USUM2/RF ( IFP  ) 

GO  TO  1060 
1019  CONTINUE 

IF(RF( IFP). LE. 1.0  )  GO  TO  1016 
GO  TO  1016 
1016  CONTINUE 

C**  S  -  R  SOMEWHERE - RESTRICTION 

ISI  -  UR(IFP) 

UR(IFP)-1.0 

UX0»0.5*DELPS(ISI ) / { RHOF+US ( 1 S I ) ) 

UX1*0.0 

DO  1014  IXI«1#ISI 

UX1  •  UX1  +  DELPS(1XI)*S(IXI)/(US(1XI)**2) 

1014  CONTINUE 

IF  C ISI.EQ.1.0R.ISI.EQ.16)  GJ  TO  93 
UPXP  »  (LS(ISI+1)-US(ISI-1) )/(2.*DS) 

URX(IFP)  ■  UXO  ♦  DS*URXP*UX1/(PF(  IFP)*PHOF) 

93  CONTINUE 

UX1  ■  UX1*AWK*UR( IFP)**2*0.5/(PF( IFP)**2*RH0F) 

UX1  «  UX1*DS*(US< ISI )**2) 

UR(IFP)  »  US(ISI) 

U(IFP)  »  UX0-UX1 
GO  TO  1060 
1016  U(IFP)  ■  0.0 
1060  CONTINUE 
C** 

C+*  SAMPLE  CASE  THRUST  COEFFICIENT 

C** 

U(IFP)  -  U(  IFP)/CTHRUST 
URX(IFP)  «  URX( IFP)/CTHRUST 

C** 

lCul  CONTINUE 
C 

C  PRINT  RESULTS 

C 

CALL  OUTPUT 

2001  CONTINUE  ,  % 

STOP  (B2) 


oooooooo 


END 

SUBROUTINE  INPUT 

COMMON/I NPUT/RHOF, S<32 ), DS»XF( 100),  RF (100  >,DELPS< 32) 

CQMM0N/NCMWAKE/US(32),N0FPS»0R(100),U(100),AWK 

COMMON/LAGZRS/  X  C 15 ) >  N 

COMMON/STRIP/  NSTRIP 

CON MON /THRUST /  C THRUST* UNORM 

READ  INPUT  —  RHGF-FLUID  0 ENS I T Y, S ( I )  «  I  STRIPS,  US-WIDTH  OF 
STRIPS#  ( XF , RF  )  “COORDS •  OF  FIR 

STRIPS#  (XF,RF)-COORDS.  OF  FIELD  POINTS,  WHERE  XF.LT.Q.O,  Ui(I)» 
INFLOW  AT  EACH  STRIP 

INPUT  THE  MOD  OR  ABSOLUTE  VALUE  OF  XF 

DIMENSION  IDENT ( 7  ) 

READ  (5,100)  ( IDENT(I), 1-1,7) 

WRITE (6, 101 ) 

WRITE( 6#  212) 

WR  ITE  (6,  100  )  (IOENT(I ), 1-1,7) 

READ (5, 102)  RHOF,DS,AWK,C THRUST, UNORM 
IF(CTHRUST.EQ.O.O)  CTHRUST-1.0 
IF(UNDRM.EQ.O.O)  UN0RM-1.0 
WRITE (6,213 ) 

WRITE (6, 110) 

WRITE! 6, 111 )  RH0F,D$,AWK 

RE  AD ( 5, 103 )  NSTRIP 

RE AD ( 5, 102 )  (S(I),I-1, NSTRIP) 

WR I TE ( 6, 1 12  ) 

WR I TE ( 6, 102 )  (S(I), 1*1, NSTRIP) 

READ (5, 102 )  (LS(I), 1-1, NSTRIP) 

WRIT  E ( 6,  113 ) 

DO  1001  IXI-1, NSTRIP 
1001  US(IXI)  •  US ( 1X1 ) /UNORM 

WRITE ( 6, 102 )  (US( I), 1-1, NSTRIP) 

READ (5, 102)  ( DELPS ( I), I-l, NSTRIP) 

WRITE ( 6, 114 ) 

WR I T  E( 6, 102 )  (DELPS( I), Iai, NSTRIP) 

RE  AD (5, 103)  NOF  PS , N 

IF(NOFPS.GT. 100)  WRITE(6,104)  NOFPS 
IF(NQFPS.GT.IOO)  STOP 
WR  IT  E ( 6, 1 15 )  NOFPS 
WRI TE ( 6, 116  ) 

DO  105  1*1, NOFPS 

105  RE  AD ( 5, 102  )  XF ( I ) , RF ( I ) , UR ( T  ) 

DO  117  1*1, NOFPS 
XFA  -  -XF(I) 

117  WRITE(6,118)  XF A, PF ( I ) , UP ( I  ) 

C 

C 

C  FORMAT  STATEMENTS  FOP  INPUT/OUTPUT 

C 

101  FORMAT ( 1H1 ) 

100  FORMAT ( 7A10 ) 

212  FORMAT!//) 

213  FORMAT!///) 

103  FORMAT (815) 

110  FORMAT ( /10X, *  INPUT  CARD  NO.  2*/) 

111  FORMAT ( 5X,*  F HOF»*F 10 .5, *  DS«*F10.5»*  AWK--F10.5/) 

112  FORMAT ( /10X, *  INPUT  CARD  NO.  3, PROPELLER  STR i PS-M 1 DPO INT S* / ) 

102  FORMAT ( 8 F10. 5  ) 

113  FORMAT ( / 10X, *  INPUT  CARD  NO.  4, INFLOW  AT  EACH  STRIP*/)  (83) 


F0RMATC/10X,* 

LIP*/) 

FORMAT ( / /10X* *  NO.  OF  FIELD  POINTS  **I3) 

FORMAT ( / 10X#  *  FIELD  POINT  COORDINATES*/) 

FORMAT ( 10X>*  XF  **F10.5,8X,*  RF-*F10.5,*  UR«*F10.5) 

FORMAT ( / 5 X» *  NOFPS«*I5>*  EXCEEDS  DIMENSIONS  OF  100*/) 

UR  IS  THE  VELOCITY  OF  STREAM  AT  RF 

RETURN 

END 

SUBROUTINE  INTGRD ( IFP ) 

COMMON/ I NTGkND/FI0>  FIX* FI2»  XIO* XI1>  XI2 
DIMENSION  FIO(960),FI1(960),FI2(960 ) 

DIMENSION  XI0(32)#Xll(32)#XI2( 32) 


INPUT  CARD  NO. 
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5,  STEADY  THRUST  L0A01N6  AT  EACH  ST  I 


CALCULATION  OF  INTEGRANDS  FOR  I1,I2,AN0  10 

COMMON/BESFUNC/XJOS* XJ0R>  XJ1R 
DIMENSION  XJ0S(960)*XJ0RC50)>XJ1R(5Q) 

COMMON/ INPUT /RHOF>S<  32) »DS#  X  F  < 100)* RF< 100), DEL  PS  (  32) 

C0MM0N/N0MWAKE/US<32)*N0FP$»UR(100),U( 100), A WK 

COMMON /LAGZRS/  Y(15),N 

COMMON/LAGZRS/  YY(15),NNN 

C0MM0N/LAG/XL(50),AL(50),NN 

COMMON/STRIP/  NSTRIP 

DIMENSION  XJ(10),XJ1<10) 

DIMENSION  Y ( 50 ) 

IF(N.EO.O)  N *15 

GET  ZEROS  OF  LAGUERRE  POLYNOMIALS 

CALL  LGZEROS(N) 

N*NN 

DO  300  I *1*  N 
Y( I)-XL(I) 

IF(N.EQ.10.0R.N.EQ.15 )  CALL  LG'ZEROS(N) 

CALCULATE  BESSEL  FUNCTIONS  JO  AND  J1 

DO  100  I  *1,  N 
DO  101  J  *  1>NST  R I P 
IJ  -  NSTR IP* ( 1-1 )  + J 
ARGO ■ Y ( I)*S(J)/XF(IFP) 

CALL  BESLC1,XJ,0,10,0,ARG0,1) 

XJ0S(1J)-XJ(1) 

CONTINUE 

ARG1*Y( I)*RF( IFP)/XF(  IFP) 

CALL  BESLC1#XJ1#0>1G,0>AR61,2) 

XJOR(I)-XJKl) 

XJ1R(I )-XJl(2) 

CONTINUE 

DO  200  J-l.NSTRIP 

DO  201  I  *1,  N 

IJ  ■  NSTRIP*( 1-1) *J 

CALCULATE  FIO 

FIOC IJ)*Y(I)*XJOS( IJ )*XJOR( 1)/(XF(1FP)**2) 

CALCULATE  FI1  AND  FI2 


FI1C  IJ)-XJ0S(IJ)*XJ1RCI)/XF(IFP) 


FI2(IJ)  •  (1.-XJ0R(I))*XJ0S(1J)/(Y(I)*RF(IFP)) 
201  CONTINUE 
200  CONTINUE 
RETURN 
END 

SUBROUTINE  IGZEROS(N) 

COMMON/ LAG ZRS /  X(15),N1 

ZEROS  OF  THE  LAGUERRE  POLYNOMIALS 

SELECT  N  -  10/  OR  N  ■  15  ONLY 


DIMENSION  XIO(IO) 
DATA ( ( X10 ( I ) , I *1, 
11.80834290174,3.4 
28.330152746764/11 

3  21.996585811981/ 
DATA  ((X15(I),I«1 

1  1.215595412071/2 

2  5.425336627414/7 
3/13.130262482176/ 

4  25.623894226729, 

5  48.026085572686) 


/  XI 5 ( 1 5  ) 

10 )«.  13 7793470540/ .7294  545  49503/ 
014336977855/5.5552496140064/ 
.8437856379,16.279257631378, 
29.920697012274) 

,15) *0.093307612017, 0.49269 1740302/ 
.269949520204,3.667622721751/ 
.565916226613/10.120228568019 
16.654407708330/20.776478699449, 
31.407519169754/ 38.530683306486/ 


1 

■4 


N1  -  N 

DO  100  1*1/15 
X(I)  «  0.0 

IF(N.EQ.15)  GO  TO  102 
IF(N.EQ.IO)  CONTINUE 
DO  101  1*1/10 

xm-xiom 

GO  TO  104 
CONTINUE 
DO  103  1*1/15 
X(I)  -  X15( I) 

CONTINUE 

RETURN 

END 

SUBROUTINE  GAUSSLG 

C0MMGN/INTGRND/FI0,FI1/F12,XI0,XU,XI2 
DIMENSION  FI01960 )/ FI 1(960) / FI2( 960) 
DIMENSION  XI0(32)/X11(32)/X 12(32) 
DIMENSION  A(50)/A10(10),A15(15 ) 

COMMON /LAG/ XL (50),AL<  50)/ NN 
COMMON/JTR IP/  NSTRIP 

GAUSS-LAGUERRE  QUADRATURE 

COMMON /LAG ZRS/  X(15),N 


DAT A( ( A10( I)/ 1*1/ 10) *.308441115765/ .401119929155/ 

1  .218068287612/0.0620874560987/ .00950151697518, 

2  .000753008388588/ .000028259233496/ .0000004249314/ 

3  .000000001839565/ .000000000000991183) 

DATA  ( ( A15 ( I )/  I *1/15 >*.2182 348 8594/ .3422101779  23/ 

1  .263027577942/ . 126425818106/ . 040206864921/ . 0o8 5638 77&C 361/ 

2  .00121243614721/ .000111674392344/ .00000645992676202/ 

3  .00000022263169071/ .00000000422743038498/ 

4  .000000000039218973/ .000000000000145652/ 

5  .000000000000000146303/ .000000000000000000016006 ) 

IF(N.E0.55)  GO  TO  500  ,  . 

N-NN  <B5> 
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IF(N.NE.10.0R.N.NE.15)  GO  TO  106 
IF(N.EQ.IO)  CONTINUE 
IFIN.EQ.15)  GO  TO  101 
DO  100  1-1*10 
A ( 1 1  *  A10(I) 

GO  TO  103 
DO  102  1-1*15 
A(I)  -  A 1 5  f I ) 

GO  TO  103 
DO  108  I-1#N 
Ad)«ALd) 

CONTINUE 

XI0#XI1»XI2  INTEGRALS 


DO  10*  J -1*  NSTR  I P 
SUMO-0.0 
SUM1-0.0 
SUM2-0.0 
DO  105  I»1*N 
IJ»NSTRIP*(I-1 )+J 
SUMO •SUMO* FI 0 ( I J) *Al I ) 
SUM1»SUM1+FI1(IJ)*A(I ) 

SUM2-SUM2+FI2 (IJ)*A(I) 

105  CONTINUE 

XIO(J)  -  SUMO 

C**  WRITE (6/220)  XIO(J) 

220  FOPMAT ( / /*  XI 0-*E 1*. 6/ / ) 

XI1(J)  »  SUM  1 
XI2(J)  -  SUM2 
10*  CONTINUE 
GU  TO  505 

500  CONTINUE 
N  ■  15 
SUMO-0.0 

DO  501  I-l>15 

SUMO  -  SUMO  ♦  FI0d)  +  Al5d) 

C**  WRITE (6* 5050)  SUMO* F 1 0( I ) * A15 ( I ) 

5050  FORMAT ( 2  X*  ♦  5UM0-*fel* .6* *  F10«*E1*.6»* 

501  CONTINUE 

X 10 ( 1 )  -SUMO 

WRITE ( 6*  220 )  XI0(1) 

RETURN 

505  CONTINUE 
RETURN 
END 

SUBROUTINE  BESL(M#J*Y*NDIMJ#NDIMY#X*NMAX) 
C 

C  CALCULATION  OF  BESSEL  FUNCTIONS 

C  M«l»  FOR  J 

C  M-2  *  FOR  Y 

C  M-3  *  FOR  J  AND  Y 

C  NDIMJ-DIKENSION  OF  JUT  LEAST  2) 

C  NOIMY-DIMENSION  OF  Y 

C  X-ARGUMENT 

C  NMAX-MAXIMUM  ORDER  PLUS  ONE 


■  FI0d)*Al5d) 
l  SUM0*FI0d)*Al5d) 

SUM0-*E1* .6* *  F10-*E1*.6»* 


XI0( 1) 


Al 5 1 1 ) -*  E 1*  •  6 ) 


J ( AT  LEAST 
Y 

PLUS  ONE 


R  £  A  L 

DIMENSION  J(ND1MJ)*YCNDIMY)#CJ0L3I7)*CJ1L3(7)*CY0L3(7)»CY1L3(7)» 
1CF0(7)»CTH0(7)#CF1(7)*CTH1(7> 

DATA  t(CJ0L3d)*I-l#7)-  1 .0* -2 .2*99997#  1 . 2656208*- . 3163b 660# 

1  .0****79C»-. 0039***0#  .00021000)* ( I CJ 1L 3 ( I ) # I- 1* 7 ) - . 50  # 


-1 ********  - 
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2-. 56249985/  . 2 109 3573, -. 0395 4289,  .00443319/ -.00031761/  .00001109  ) 
3/  ((CY0L3(I), 1-1,7)-  .36746691,  .  605  59  3t>b,  - .  74  3  5038  4,  .25300117/ 
4-. 04261214/  .0042 7916/ -.000  248 46)/  (  (CY1L3(I), I -1, 7 ) — . 63061980, 

5  .2 2120910/  2. 166 2709/ -1.31 64 327/  .31239510/-. 04009760/  .0027673) 
DATA  ( (CFO (I), 1-1,7)-  .7976 8456 /-. 00000077/-. 005 5274 0/ - .00009512/ 

1  .00137237/-. 00072805/  .00014476)/  C CCTHO ( I >, 1 -1/ 7) — . 7853931b/ 

2-. 04166397/-. 00003954/  . 00262513/ -.000541 25/-.00029333, .00013558 ), 
3( (CF1( I)/ 1-1/7)-  .79788456/  .00000156/  .01659667,  .00017105/ 

4-. 00249511/  .00113653/-. 00020033  )/  ((CThl(I), 1-1,7)*- 2.3561 9<*  49, 

5  .12499612/  .00005650/ -.00637879/  .00074348/  .00079824, 
fc-. 00029166) 

IF ( X. EQ.O . )  60  TO  9 
PI-  3.1415926 
Xl-ABS(X) 

S-X/Xl 
T-  X*X/9 . 0 

IF(M.NE.l)  FLPI-2.*A10G(.5*X1)/PI 
IF ( XI . GT • 3 )  GO  TO  1 
J(l)»  P0LY(6,CJ0L3,T,7) 

IF  <M  .NE.  1)  Y(1)-P0LY(6,CY0L3,T,7)+J(1)*FLPI 
IF  ( NMAX  .EQ.  1)  RETURN 

IF  ( M  .EO.  1  .AND.  NMAX  .GT.  2)  GQ  TO  5 
J(2)-X+P0LY(6,CJ1L3/T,7) 

IF  (M  .EQ.  1)  RETURN 
Y ( 2 ) -  FLPI*J(2)+PQLY(6,CY1L3/T,7)/X 
IF  (NMAX  .EQ.  2)  RETURN 
I F ( M  .EQ.  2)  GO  TO  6 

5  MM-NMAX+5 
FM-MM 

J ( MM+2  ) • 1 .0 

J ( MM+1 )  *  SORT  ( FM* ( FM  +  1  •  )  )*(  (  FM+l  .  )  /FM )  *♦  ( MM+1 )  *2  •  /  (  X*2 . 71828  ) 
A-J(l) 

MM-MM+2 

CALL  BESRECd,  J/NOIMJ/X/MM) 

FAC-  A/J(l) 

DO  7  I-l/NMAX 
7  J(I)-FAC*J(I) 

IF  (M  .EQ.  1)  RETURN 

6  CALL  BESREC(2, Y/NDIMY/ X/NMAX  ) 

RETURN 

1  T-3./X1 

SQX-1 . /SORT ( XI ) 

F-P0LY(6/CF0/T/7) 

TH-P0LY(fe/CTH0/T/7)+Xl 
IF  (M  .NE.  2)  J(1)-F*C0S(TH)*SQX 
IF  (  M  .NE.  1)  Y(1)-F*SIN(TH)*SQX 
IF  (NMAX  .EQ.  1)  RETURN 
FNMAX-NMAX-1 

IF  (M  .EQ.  1  .AND.  NMAX  .GT.  2  .AND.  FNMAX  .GT.XDGO  TO  5 
F-  POLY (6/CF1, T/7) 

TH-  P0LY(6/CTH1/T,7)+Xi 

IF ( M  .NE.  2  .AND.  FNMAX  .LE.  XI ) J ( 2 ) -F *C QS ( TH ) * SQ X* S 
IF  (M  .NE.l)  Y(2)-F*SIN(TH)*SQX 
IF  (NMAX  .EQ.  2)  RETURN 
IF  (M  .EQ.  2)  GO  TO  6 
IF  (FNMAX  .GT.  X1)GC  TO  5 
CALL  BESREC(2, J/NDIMJ, X, NMAX ) 

IF  (M  .EQ.  3)  GO  TO  6 
RETURN 

9  DO  10  I-l/NMAX 
10  J(I)-0. 

J ( 1 )-l. 


(B7) 


) 

) 

) 


) 


c 


) 

) 

) 

J 

J 


RETURN  9 

END 

SUBROUTINE  BE SREC (M, U, ND IM 4, X, NM AX ) 

RECURRENCE  FORMULAS  FOR  BESSEL  FUNCTIONS 

M*l,  FOR  BACKWARD  J  AND  Y 

N«2,  FOR  FORWARD  J  AND  Y 

M-3,  FOR  BACKWARD  I 

M-A,  FOR  FORWARD  K 

W(I)»  INDICATED  BESSEL  FUNCTION  OF  ORDER  1-1  AND  ARGUMENT  X 
NMAX*  LARGEST  ORDER  PLUS  ONE 


DIMENSION  W(NDIMW) 

F1(FN, U,V)*  2 »  *FN  *U/ X— V 

F2(FN,U,V>-  2 • *FN  *U/ X  +  V 

IF  (M  .EQ.  1  .OR.  M  .EQ.  3)  GO  TO  1 

IF  ( M  .EQ.  A)  GO  TO  6 

DO  3  1-3, NMAX 

FI*  1-2 

3  W<I)«  Fl(FI,W(I-l),W(I-2>) 

RETURN 

1  IF  (M  .EQ.  3)  GO  TO  A 
DO  2  1*3, NMAX 

J-  NMAX-I+1 
F  J  *  J 

2  W  ( J  )  *  F1CFJ,W(J+1  ),W(J+2)  ) 

RETURN 

A  DO  5  I*3,NMAX 
J-NMAX-I+1 
F  J  *  J 

5  W C J ) *  F2(FJ,W(J+1  )  ,W ( J  +  2 ) ) 

RETURN 

6  DO  7  1*3, NMAX 
FI-I-2 

7  W(I)*F2(FI,W( 1-1) » W( 1-2 ) ) 

RETURN 

END 

SUBROUTINE  MODBE S ( M, 1 , K» ND IM If NDINK, X,NM A X  ) 

COMPUTATION  OF  MODIFIED  BESSEL  FUNCTIONS 

M*l,  FOR  I 

M«2*  FOR  K 

M*3,  FOR  I  AND  K 

X* ARGUMENT 

NMAX-MAXIMUM  ORDER  PLUS  ONE 


REAL  I,K 

DIMENSION  I(NDIMI),  K(NDIMK),  CICL  (7),  CI1L  (7)/  CKOL  (7), 

1CK1L  ( 7 ) »  CIOG  (9),  CI1G  (9),  CKOG  (7),  CKIG  (7) 

DATA  ( (CIOL( I), I*l,7)*l.,3.5156229, 3.0899A2A,1 ,20o7A92,  .2659732, 

1  .0360766, .00A5813), ( ( C 10G ( I ) , I • 1, 9  )  * .  3989A226 , .01326592, 

2  .00225319, -.00157565, .00916261, -.02057706, .02635537, -.016A7633, 

3  .00392377), ( ( C 1 1L ( I ) ,  1*1, 7 ) * . 5, .6789059A, .51 A98669, . 1 50B A9 3A, 

A  .02658  733, .00301532,  .000  32A11), <(CI1G(I), 1*1,9)*. 3989A2 28, 

5  -.0398802 A, -.00362018, .00163601,-. 01031 555, .022 82967, -.028 953 12, 

6  .0178  765 A, -. 00 A  20059 ) 

DATA  (  (CKOL  (I),  I*  1,7)— .57721566,  .A2278A20,  .23069756,  .03*86590, 

1  .00262698, .00010750, .000007A ),( (CKOG ( I), 1*1, 7) *1.2533 1A1A, 

2  -.07832358, .02189568,-. 01062 AA6, .00587872, -. 002 5 15 A, . 0005 320tt ) , 

3  ( (CKIL(I), 1*1,7)*  1., .15AA31AA, -.67278 579, -.18 156897, -.019 19 A02, 
A  - . 00110A0A,- .0000A68 6 ) , ( ( CKIG ( 1 ) , I *1, 7 ) *1 • 25 33 1 A1 A, • 23A98619, 

5  -.0365562, . 0150A268,-. 00780353, . 0032561 A, -.000662 A5 ) 

GCOEF  (  FN,  X  )  *  SQRT(FNAFN  +  X*X)+FN*AL0G(X/(FN4SQf«T(FN*FN  +  X*X)  )  ) 

IF  (X  .LE.  0.  .ANO.  M  .NE.  1)  60,10  200 


J 
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t 

s 

S 

) 

) 

) 

) 

) 

) 

I 

) 

) 

) 

) 

» 

) 

) 

I 

) 


IF  (X  .EQ.  0.)  60  TO  99 
XI-  ABS(X) 

T2-  .25*X*X 

TR-2./X 

Tl-X/3.75 

T1S-  X*X/14.0625 

T1R-  3.75/X1 

XSR-  I » / SORT ( XI ) 

DEC-  EXP(Xl) 

DECK-  EXP (  -XI) 

FIN-  ALOG ( • 5  *X 1 ) 

S-  Xl/X 

IF  ( M  .EG.  2)60  TO  100 
C 

C  MODIFIED  BESSEL  (FIRST  KIND) 

C 

IF  (XI. GT.  3.75)  GO  TO  10 
I(l)«  PQLY(6»CI0L*T1S>7) 

IF  (M  .EO.  1  .AND.  NMAX  .EQ.  1)  RETURN 

IF  ( M  .EQ.  3  .AND.  NMAX  .EO.  1)  GO  TO  100 

IF  (NMAX  .NE.  2)  GO  TO  1 

1(2)-  P0LY(6,CI1L>T1S»7)*X 

IF  (M  .NE.  1)  GO  TO  100 

RETURN 

1  A-I ( 1 ) 

MM-  NMAX*5 
FM-  MM 

I (MM  +  1 ) -  l./SQRT(6. 28316  *FM)/((1.+X*X/(FM*FM) ) **0 • 25 ) * 

1EXP(GC0EF(FM*X1))*S 
FM-FM+1 

I (MM+2 ) -  l./SQRT(6. 28318  *FM ) / ( ( 1 . *X*X/ ( F R *F K ) ) **0 . 25 ) * 
1  EXP ( GCOE f ( FM* XI)  ) 

MM-  MM+2 

CALL  BESREC(3,I>NDIMI,X,MM) 

FAC-  A/I ( 1) 

DO  2  L-l.NMAX 

2  I ( L  )  -  F AC* I ( L  ) 

IF  (M  .EQ.  3)  GO  TO  100 
RETURN 

10  1(1)*  PGLY(  8,CI0G,Tllo9)*XSR*DEC 

IF  (M  .EQ.  1  .AND.  NMAX  .EQ.  1)  RETURN 
IF  (M  .EQ.  3  .AND.  NMAX  .EQ.  1)  GO  TO  100 
IF  (NMAX  .NE.  2)  GO  TO  1 
1(2)-  P0LY(8,CI1G,T1R#9)*XSR*0EC*S 
IF  (M  .NE.  1)  GO  TO  100 
RETURN 
C 

C  MODIFIED  BESSEL  (SECOND  KIND) 

C 

100  IF  (X  .GT.  2.)  GO  TO  110 

IF  (M  .EQ.  2)  I ( 1 ) «  PULY(6»CiOL»TlS»7) 

K ( 1 ) «  -FLN*I(l)*P0LY(6/CK0L>T2/7) 

IF  (NMAX  .EQ.  1)  RETURN 

IF  (M  .EQ.  2)  1(2)-  POLY  ( 6, C IiL# T1 S> 7 ) ♦ X 
K ( 2 ) -  <X*FLN*I(2)+PGLY(6,CK1L#T2#7) )/X 
IF  (NMAX  .NE.  2)  GO  TO  101 
RETURN 

101  CALL  BESREC(4*K#NDIMK,X#NMAX) 

RETURN 

110  Ml)-  P0LY(6,CK0G,TM7)*XSk*DECK 
IF  (NMAX  .EQ.  1)  RETURN 
K (2  )  -  P0LY(6»CK1G,TR>7)*XSR*DECK  ,  ^ 
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IF  ( NMAX  .NE.  2)  60  TO  101 
RETURN 

99  00  98  J«1#NMAX 
98  I(J)>0. 

I(1I>1. 

RETURN 

200  PRINT  97#  X 

97  FORMAT  (10X26HK  NOT  DEFINED  FOR  X  .LE.  0# 2 X2HX «E 16 . d  ) 

CALL  EXIT 
END 

FUNCTION  POLY(N,A#X#NDIMA) 

DIMENSION  A(NDIMA) 

POL  Y*A ( N+l ) 

DO  1  I  *1#  N 
J-N+l-I 

1  POLY-Ai J)*X*POLY 
RETURN 
END 

SUBROUTINE  OUTPUT 
C 

C  WRITE  OUT  THE  COMPUTED  VELOCITIES  AT  THE  GIVEN  FIELD  POINTS 

C 

C0MM0N/INPUT/RH0F#S(32)#DS#XF< 100 ) # RF < 100  ) # DE LPS ( 32 ) 

COMMON/ NOMWAKE /US ( 32 ) #  NQF  PS# UR ( 100 ) #  U( 100 ) »  AWK 
COMMON/DUOR/  URX(IOO) 

WR ITE ( 6# 100  ) 

100  FORMAT ( 1H1 ) 

WR ITE ( 6# 101 ) 

101  FORMAT t/10X#*  VELOCITY  AT  GIVEN  FIELO  POINTS*/) 

WR ITE ( 6# 102  ) 

102  FORMAT ( / 5X#  *  ABS(XF)  RF  UR  U  */) 

WRITE (6# 105)  (XFC I ># RF ( I ) #UR( I )# U < I ) # 1-1, NCFP S  ) 

105  F0RMAT(5X#F10.5#2X#F10.5#2XF10.5#5X#E14.6) 

WR I TE ( 6# 100  ) 

WRITE ( 6# 110) 

WRITE ( 6# 102  ) 

WRITE (6# 105)  (XF( I ) # R F ( I ) # UR ( I ) # URX ( I ) # I *1# NOF P S ) 

110  FORMAT ( / 10X#  *  VELOCITY  FIELD  ON  DISK  USING  FXACT  DUDR*/ ) 
RETURN 
END 

SUBROUTINE  L AGUER < NN# X, A# AL F# B#C> E PS#C SX# CS A# T S X# TS A ) 
DIMENSION  X(50)#A(50)#B(50)#C(50) 

FN  -  NN 
CSX  •  0.0 
CS A  -  0.0 

CC  -  GAMMAt  ALF+1 • ) 

TSX  -  FN*( FN+ALF ) 

TSA»CC 
DO  1  J«2»  NN 

1  CC  ■  CC*C(J) 

DO  7  I *1#  NN 
IF(I-l)  6#  2#  3 

2  XT  •  (1.*ALF)*13.*.92*ALF  )/ (1.+2. 4*FN+l.d*ALF ) 

GO  TO  6 

3  IFtI-2>6#4#5 

4  XT  ■  XT*(15.*6.25*ALF ) / ( l.*.9*ALF  +  2»  5*FN ) 

GO  TO  6 

5  FI  ■  1-2 

R1»(1.*2.55*FI)/<1.9*FI) 

R2-1.26*FI*ALF/<1.*3.5*FI) 

RATIO* (R1*R2 )/(l.*.3*ALF) 

XT  •  XT  ♦  RATIO*  C  XT-A  ( 1-2  )  )  (BIO) 
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6  CALL  LG ROOT (XT#NN#ALF#DPN#PN1#B#C#EPS) 

X(I)  -  XT 

All)  •  CC/DPN/PN1 
C**  PUNCH  20#  ALF#  NN# 1#  XT#  A  1 1 ) 

CSX  »  CSX  ♦  XT 

7  CSA  »  CSA  ♦  All) 

C*+  PUNCH  20,  ALF#NN, I#CSX,CSA# TSX, TSA 

20  F0PMATIF6.2»2I3#2(1X#E14.8),2X#2<1X#E14.8  )  ) 

RETURN 

END 

SUBROUTINE  LG ROOT (X#NN#ALF#DPN#PN1#B#C#EPS) 

DIMENSION  B( 50 ) #C I  SO ) 

ITER  «  0 

1  ITER  »  ITER  +1 

CALL  LGRECR(P#DP#PN1#X#NN#ALF#B#C) 

D  -  P/DP 
X  ■  X-D 

IF(ABS(D/X)-EPS)  3#  3#  2 

2  IF ( ITER-10)  1#  3#  3 

3  DPN»DP 
RETURN 
END 

SUBROUTINE  LGR ECR I PN# DPN# PN1# X# NN# AL F # B# C  ) 

DIMENSION  B ( 50  )  # C I  50 ) 

Pl-1. 

P  •  X-ALF-1. 

DP1-0. 

DP  ■  1. 

DO  1  J«2#NN 

0  «  IX-BIJ) )*P  -  CIJ  )*PI 
DQ  ■  I X— B I J ) ) *DP  ♦  P  -  CIJ)*0P1 
PI  «  P 
P  «  0 
DP1-DP 
1  DP  -  DQ 
PN  «  P 
DPN  ■  DP 
PN1-P1 
RETURN 
END 

FUNCTION  GAMMA! X) 

GAM(Y)  -  (((((((. 035868343*Y-.1935278ie)*Y+.482199394)*Y 

1  .756704078 ) *Y+.918206857 )  *Y-. 897056937 )* Y+ . 96 8205891)* Y 

2  -.577191652 ) *Y  ♦  1.0 
Z  ■  X 

IF ( Z )  1# 1#4 

1  GAMMA  >  0.0 
C**  PUNCH  2#  Z 

2  FORMAT (2X#19HARG  ERROR  FOR  GAMMA  #  E16.8) 

WRITE(6#2)  Z 

GO  TO  14 

4  IFIZ-70.)  6# 1# 1 

6  IF ( 2-1 . )  8#  7# 9 

7  GAMMA  >  1.0 
GO  TO  14 

8  GAMMA  *  GAM ( Z ) / 2 
GO  TO  14 

9  ZA-1.0 

10  Z  -  Z-l. 

IFIZ-l.)  13# 11# 12 

11  GAMMA  •  ZA 

GO  TO  14  (B It) 


****** 


.1  srff* 


LESS  THEN 


EQUAL 


12  ZA  •  Z  A*  2 
GO  TO  10 

13  GAMMA  -  ZA*GAM(Z) 

14  RETURN 
END 

SUBROUTINE  LEG2HQ(X,N,Q,Q1,IC) 

A  FORTRAN  SUBROUTINE  TO  EVALUATE  LEGENDRE  FUNCTIONS  OF 

KIND  AND  THEIR  DERIVATIVES 

X-ARGUEMENT 

N-3/2-MAXIMUN  ORDER  REOUIRED  (N  LESS  THEN  OR  EQUAL  TO 

Q-LEGENDRE  FUNCTION  OF  SECOND  KIND 

Ql-DERIVATIVE  OF  Q 

IC»1  FOR  JUST  Q 

IC-2  FOR  0  AND  Ql 

DIMENSION  Q<125),Q1(125>,R(125> 

11  M-N+19 

6  FM»M 

IF(X-1.0C14)  7,8,8 

7  QCM+1)-RLEG2(X,M+1,19) 

Q(M>»RLEG2JX,M,19) 

GO  TO  9 

8  0<M)-1. 

Q(M+1)-(FM*X-SQRT(FM*FM*(X+1»)*< X“1 •)+«25))/(FM+«5) 

9  DO  1  I J-2»M 
I-M-IJ+2 

FI *1-1 

F-1.0/ (FI-0.5) 

Q(I-1)-2»*FI+X*F*Q(I)-(F1+0»5)*F*Q( 1*1) 
IF(Q(I-1)-1.E+30)1,1*10 
10  N-N-10 
GO  TO  11 

1  CONTINUE 
S«ALEG1(X,1> 

S-S/Q( I) 

DO  2  I«1,N 

2  Q(I)-S*0(I) 

IF(IC-l)  4,4,3 

3  S-l./(  (XU.)MX-l.)  ) 

<J1U>— .5*SMX*Q(l)-0(2)) 

DO  5  1-2, N 

FI- 1*1 

5  Q1(I)-(FI— •5)*S*(X*Q(l)-Q(I-l) ) 

4  RETURN 
END 

FUNCTION  RLEG2 ( Z, N,M ) 

QAR6<U,FN,W,W2  >»( ( W/ < W+U > ) **FN-1 . ) / SORT < U* ( U*W2 ) ) 
W-Z+SQRT( <Z+1. )*(Z-1. ) ) 

W2-2.MW-Z) 

FN-N 

FN-FN-.5 

RLEG2-0. 

U-.001 
DO  1  1-1,19 

RLEG2-RLEG2  +  UAR6(U,FN,W,  W2) 

1  U-U+.001 
A-QARG(»02,FN,W, W2) 

RLEG2-(.5>*A+RLEG2)*.001*A*.005 

U-.03 

00  2  1-3,  M 

RLEG2-RLEG2*.01*QARG(U,FN,W,W2) 

2  U-U+.01 

RLEG2-RLEG24. 005*QARG( U, FN, W, W2 ) 
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A-SQRTC (Z-l. >*CZ+1.  )  ) 

B— AlOGC  A) 

C-ALOG(U*A*SQRT(U*(U*2.*A) ) ) 

D-RLEG2+B+C 

RLEG2«D*Vi**(-FN) 

RETURN 

END 

FUNCTION  ALEG1CX/N) 

DIMENSION  Z(20)/A(10)/B(10) 

PI-3.1415926 
I F  {  X—l  •  09 )  3/4/4 

4  DO  1  1-1/20 
FI-D2-1 

1  2(I)-C0S(FI*PI*.025) 

ALEG1-0.0 

FN-N 

DO  2  1-1/20 
A1«1.0-Z(I)*Z(I) 

A2- X-Z ( I ) 

2  ALEG1-ALEG1+(A1/A2)**{N-1>/SQRT< A2) 
ALEGl»ALEGl/(2.0**( FN-O. 5) )*0.05+PI 
RETURN 

3  Y-X-1.0 
AL-ALOGCY) 

A(1)-2.5*AL0G(2.) 

IF(N-l)  5/5/6 

6  NZ-N-1 

DO  7  1-1/ NZ 
F I  *  I 

7  A(1)»A(1)~2./(2.0*FI-1*) 

5  BCD  — .5 
FN-N-1 
G-FN*FN-.25 
DO  9  1-2/10 
FI-I-2 
H-I-l 

H2- • 5/ ( H*H) 

H3-F I*H 

Atl  >«(AU-1>MG-H3)-BC  I-1)*(2.*G/H+1.  )  )*H2 
9  BCI)»BCI-l)+( G-H3 ) *H2 
ALEG1-AC10)+AL+B( 10) 

DO  10  1-1/9 
J-10-1 

10  ALEG1-AC  J)+AL*BC  JJ+ALEGDY 
RETURN 
END 


NUMBER  OF  CASES  • 
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